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Abstract 


In  the  following  document  a  summary  is  given  of  research  carried  out  under 
award  number  FA9550-07-1-0195  from  the  Air  Force  Office  of  Scientific  Re¬ 
search.  The  objective  of  the  proposed  research  was  to  develop  and  implement 
high-order  numerical  algorithms  for  the  simulation  of  steady  and  unsteady  com¬ 
pressible  viscous  flows  with  shocks.  Notable  achievements  have  been  made  in 
three  areas,  namely  algorithm  analysis  and  development,  code  development, 
and  code  utilization  (to  investigate  various  flow  problems).  With  regards  to 
algorithm  analysis  and  development,  it  has  been  proved  for  one-dinrensional 
linear  advection  that  the  spectral  difference  method  is  stable  for  all  orders  of 
accuracy  in  a  norm  of  Sobolev  type  (provided  that  the  interior  flux  collocation 
points  are  placed  at  the  zeros  of  the  corresponding  Legendre  polynomials).  Also, 
a  new  range  of  energy  stable  high-order  methods  based  on  the  so  called  flux 
reconstruction  approach  have  been  identified.  With  regards  to  code  develop¬ 
ment,  two-dimensional  and  three-dimensional  compressible  viscous  flow  solvers 
based  on  the  spectral  difference  method  have  been  written.  The  solvers  can  run 
on  meshes  containing  straight-sided  and  curved-sided  quadrilateral  and  hexa- 
hedral  elements.  Within  the  aforementioned  codes  a  range  of  shock  capturing, 
automatic  mesh  refinement,  mesh  deformation,  and  convergence  acceleration 
algorithms  have  been  implemented  and  tested.  With  regards  to  code  utiliza¬ 
tion,  viscous  compressible  flow  over  various  two-dimensional  configurations  has 
been  investigated.  These  configurations  include  pairs  of  cylinders  (both  sta¬ 
tionary  and  rotating),  plunging  and  pitching  airfoils,  and  a  deforming  beam  in 
the  wake  of  a  cylinder.  In  three-dimensions,  turbulent  channel  flow  and  turbu¬ 
lent  flow  over  an  airfoil  have  been  investigated,  and  preliminary  simulations  of 
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viscous  compressible  flow  over  a  flapping  wing  have  been  undertaken. 


Keywords:  High- Order  Methods,  Spectral  Difference  Methods,  Shock  Capturing, 
Deforming  Meshes,  Time  Integration  Schemes 
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1  Introduction 

High-order  numerical  methods  potentially  offer  better  accuracy  than  low-order  schemes 
for  a  comparable  computational  cost.  However,  existing  high-order  methods  are  gen¬ 
erally  less  robust  and  more  complex  to  implement  than  their  low-order  counterparts. 
These  issues,  in  conjunction  with  difficulties  generating  high-order  meshes,  have  pre¬ 
vented  the  wide-spread  adoption  of  high-order  techniques  in  either  academia  (where 
the  use  of  low-order  schemes  remains  widespread)  or  in  industry  (where  the  use  of 
low-order  schemes  is  ubiquitous). 
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The  most  mature  and  widely  used  high-order  methods  (at  least  for  unstructured 
grids)  are  based  on  a  class  of  schemes  developed  in  1973  by  Reed  and  Hill  [1]  to  solve 
the  neutron  transport  equation.  Such  schemes  have  become  known  as  discontinuous 
Galerkin  (DG)  methods,  and  numerous  variants  have  been  developed  for  solving  the 
weak  form  of  both  hyperbolic  [2]  and  elliptic  systems  [3].  The  basic  principle  of 
DG  schemes  is  to  decompose  the  approximate  numerical  solution  both  spatially,  by 
tessellating  a  given  computational  domain  with  separate  elements,  and  also  spectrally, 
via  a  summation  of  piecewise  discontinuous  polynomial  basis  functions  within  each 
element.  A  particularly  simple  and  efficient  range  of  DG  schemes  utilize  high-order 
Lagrange  polynomial  basis  functions  inside  each  element,  defined  by  solution  values 
at  a  set  of  distinct  nodal  points.  Such  schemes  have  become  known  as  nodal  DG 
methods,  an  exposition  of  which  can  be  found  in  the  recent  textbook  by  Hesthaven 
and  Warburton  [4],  as  well  as  in  various  articles  by  the  same  authors  [5]  [6].  Similar 
to  nodal  DG  methods  are  spectral  difference  (SD)  methods,  (although  unlike  nodal 
DG  methods,  SD  methods  are  based  on  the  governing  system  in  its  differential  form). 
The  foundation  for  such  schemes  was  first  put  forward  by  Kopriva  and  Kolias  [7]  in 
1996  under  the  name  of  “staggered  grid  Chebyshev  multidomain”  methods.  However, 
several  years  later  in  2006  Liu,  Wang  and  Vinokular  [8]  presented  a  more  general 
formulation  for  both  triangular  and  quadrilateral  elements,  which  they  termed  the 
SD  method  (a  name  which  as  been  retained  to  the  present). 

In  the  study  presented  here,  a  broad  effort  to  investigate,  implement,  test  and  ul¬ 
timately  improve  the  SD  method  was  undertaken,  with  particular  emphasis  placed 
on  issues  such  as  shock  capturing  and  efficient  time  integration,  which  have  hitherto 
inhibited  the  widespread  adoption  of  high-order  methods  amongst  a  wider  scientific 
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community. 


2  Accomplishments 

2.1  Overview 

With  regards  to  algorithm  analysis  and  development,  the  main  accomplishments  of 
the  research  can  be  summarized  as  follows: 

•  A  proof  demonstrating  that  particular  SD  schemes  are  stable  (for  linear  advec- 
tion)  [9]. 

•  The  development  of  a  new  range  of  energy  stable  flux  reconstruction  (FR) 
schemes  [10]. 

With  regards  to  code  development,  the  main  accomplishments  of  the  research  can  be 
summarized  as  follows: 

•  The  development  of  a  two-dimensional  (2D)  viscous  compressible  SD  flow  solver. 

•  The  development  of  a  three-dimensional  (3D)  viscous  compressible  SD  flow 
solver. 

•  The  implementation  of  shock  capturing  algorithms  suitable  for  use  with  the  SD 
method  [11]  [12]. 
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•  The  implementation  of  automatic  mesh  refinement  algorithms  suitable  for  use 
with  the  SD  method  [12]. 

•  The  implementation  of  mesh  deformation  algorithms  suitable  for  use  with  the 
SD  method  [13]. 

•  The  implementation  of  convergence  acceleration  schemes  suitable  for  use  with 
the  SD  method  [14]. 

With  regards  to  code  utilization,  the  main  accomplishments  of  the  research  can  be 
summarized  as  follows: 

•  Studies  of  2D  viscous  compressible  flow  over  pairs  of  cylinders  (both  stationary 
and  rotating)  [15]  [16] . 

•  Studies  of  2D  viscous  compressible  flow  over  pitching  and  plunging  airfoils  [13]. 

•  Studies  of  2D  viscous  compressible  flow  over  a  deforming  beam  in  a  cylinder 
wake. 

•  Studies  of  3D  turbulent  channel  flow  [17]. 

•  Studies  of  3D  turbulent  flow  over  airfoils. 

•  Preliminary  studies  of  3D  flow  over  a  flapping  wing. 

Journal  articles  resulting  from  the  research  include: 
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•  P.  E.  Vincent,  P.  Castonguay.  A.  Jameson.  A  New  Class  of  High-Order  En¬ 
ergy  Stable  Flux  Reconstruction  Schemes.  Submitted  to  Journal  of  Scientific 
Computing.  May  2010. 

•  G.  May,  F.  Iacono,  A.  Jameson.  A  Hybrid  Multilevel  Method  for  High-Order 
Discretization  of  the  Euler  Equations  on  Unstructured  Meshes,  Journal  of  Com¬ 
putational  Physics.  May  2010. 

•  A.  Jameson.  A  Proof  of  the  Stability  of  the  Spectral  Difference  Method  for  All 
Orders  of  Accuracy,  Journal  of  Scientific  Computing.  January  2010. 

•  C.  Liang,  S.  Premasuthan,  A.  Jameson.  High-order  Accurate  Simulation  of 
Low-Mach  Laminar  Flow  Past  Two  Side-by-Side  Cylinders  Using  Spectral  Dif¬ 
ference  Method.  Journal  of  Computers  and  Structures.  February  2009. 

•  C.  Liang,  A.  Jameson,  Z.  Wang.  Spectral  Difference  Method  for  Compressible 
Flow  on  Unstructured  Grids  with  Mixed  Elements,  Journal  of  Computational 
Physics.  January  2009. 

Full  length  conference  papers  resulting  from  the  research  include: 

•  P.  Castonguay,  A.  Jameson.  Simulation  of  Transitional  Flow  over  Airfoils  using 
the  Spectral  Difference  Method.  To  be  presented  at  the  AIAA  Computational 
Fluid  Dynamics  Meeting.  Chicago,  Illinois.  June  2010. 

•  K.  Ou,  A.  Jameson.  A  High-Order  Spectral  Difference  Method  for  Fluid- 
Structure  Interaction  on  Dynamic  Deforming  Meshes.  To  be  presented  at  the 
AIAA  Computational  Fluid  Dynamics  Meeting.  Chicago,  Illinois.  June  2010. 


•  Y.  Li,  S.  Premasuthan,  A.  Jameson.  Comparison  of  h  and  p  Adaptations  for 
Spectral  Difference  Methods.  To  be  presented  at  the  AIAA  Computational 
Fluid  Dynamics  Meeting.  Chicago,  Illinois.  June  2010. 

•  S.  Premasuthan,  C.  Liang,  A.  Jameson.  Computation  of  Flows  with  Shocks 
Using  Spectral  Difference  Scheme  with  Artificial  Viscosity,  AIAA  Aerospace 
Sciences  Meeting.  Orlando,  Florida.  January  2010. 

•  K.  Ou,  C.  Liang,  A.  Jameson.  A  High-Order  Spectral  Difference  Method  for 
the  Navier-Stokes  Equations  on  Unstructured  Moving  Deformable  Grids,  AIAA 
Aerospace  Sciences  Meeting.  Orlando,  Florida.  January  2010. 

•  C.  Liang,  S.  Premasuthan,  A.  Jameson,  Z.  Wang.  Large  Eddy  Simulation  of 
Compressible  Turbulent  Channel  Flow  with  Spectral  Difference  Method.  AIAA 
Aerospace  Sciences  Meeting.  Orlando,  Florida.  January  2009. 

•  S.  Premasuthan,  C.  Liang,  A.  Jameson,  Z.  Wang.  p-Multigrid  Spectral  Differ¬ 
ence  Method  For  Viscous  Compressible  Flow  Using  2D  Quadrilateral  Meshes. 
AIAA  Aerospace  Sciences  Meeting.  Orlando,  Florida.  January  2009. 

•  S.  Premasuthan,  C.  Liang,  A.  Jameson.  A  Spectral  Difference  Method  for  Vis¬ 
cous  Compressible  Flows  With  Shocks.  AIAA  Computational  Fluid  Dynamics 
Meeting.  San  Antonio,  Texas.  June  2009. 

•  K.  Ou,  C.  Liang,  S.  Premasuthan,  A.  Jameson.  High-Order  Spectral  Difference 
Simulation  of  Laminar  Compressible  Flow  Over  Two  Counter-Rotating  Cylin¬ 
ders.  AIAA  Computational  Fluid  Dynamics  Meeting.  San  Antonio,  Texas. 
June  2009. 
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PhD  dissertations  resulting  from  the  research  include: 


•  S.  Premasuthan.  Towards  an  Efficient  and  Robust  High  Order  Accurate  Flow 
Solver  for  Viscous  Compressible  Flow.  Ph.D.  Dissertation,  Stanford  University, 
March  2010. 

Awards  received  during  the  research  period  include: 

•  A.  Jameson.  Elmer  A.  Sperry  Award. 

2.2  Algorithm  Analysis  and  Development 

2.2.1  Proof  of  Stability  of  the  Spectral  Difference  Method 

It  has  been  demonstrated  that  for  the  case  of  one- dimensional  (ID)  linear  advection, 
the  SD  method  is  stable  for  all  orders  of  accuracy  in  a  norm  of  Sobolev  type  (provided 
that  the  interior  flux  collocation  points  are  placed  at  the  zeros  of  the  corresponding 
Legendre  polynomials).  The  proof  is  based  on  an  energy  method.  For  solution  poly¬ 
nomials  of  degree  k  (which  result  in  a  scheme  of  order  k+1),  stability  is  demonstrated 
with  a  norm  of  the  form 


u 


I 


(u2  +  (32kcu(k)2) 


dx, 


(2.1) 


where  [3  is  a  piecewise  constant  scaling  factor  and  c  is  a  k  dependent  coefficient.  For 
further  details  see  Jameson  [9]. 
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2.2.2  Development  of  Energy  Stable  Flux  Reconstruction  Schemes 

The  FR  approach  [18]  to  high-order  methods  is  robust,  efficient,  simple  to  imple¬ 
ment,  and  allows  various  high-order  schemes,  such  as  the  nodal  DG  method  and  the 
SD  method,  to  be  cast  within  a  single  unifying  framework.  We  have  identified  a  new 
class  of  one-dimensional  energy  stable  FR  schemes.  The  energy  stable  schemes  are 
parameterized  by  a  single  scalar  quantity,  which  if  chosen  judiciously  leads  to  the 
recovery  of  various  well  known  high-order  methods  (including  the  nodal  DG  method 
and  a  particular  SD  method),  as  well  as  one  other  FR  scheme  that  was  previously 
found  to  be  stable  by  Huynh  [18].  The  analysis  offers  significant  insight  into  why 
certain  FR  schemes  are  stable,  whereas  others  are  not.  Also,  from  a  practical  stand¬ 
point,  the  analysis  provides  a  simple  prescription  for  implementing  an  infinite  range 
of  energy  stable  high-order  methods  via  the  particularly  intuitive  FR  approach.  We 
are  currently  working  to  extend  the  formulation  to  simplex  elements.  For  further 
details  see  Vincent,  Castonguay  and  Jameson  [10]. 

2.3  Code  Development 

2.3.1  Two-Dimensional  Compressible  Viscous  Flow  Solver 

A  2D  viscous  compressible  SD  flow  solver  has  been  written  in  FORTRAN.  The  solver, 
which  is  based  on  the  approach  presented  by  Wang  [19],  can  run  on  meshes  contain¬ 
ing  both  straight-sided  and  curved-sided  quadrilateral  elements.  An  methodology 
for  applying  the  code  to  unstructured  mixed  meshes  containing  both  triangular  and 
quadrangular  elements  has  also  been  developed  [20]. 
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2.3.2  Three-Dimensional  Compressible  Viscous  Flow  Solver 

A  3D  viscous  compressible  SD  flow  solver  has  been  written  in  FORTRAN.  The  solver, 
which  is  based  on  the  approach  presented  by  Wang  [19],  can  run  on  meshes  containing 
both  straight-sided  and  curved-sided  hexahedral  elements.  To  facilitate  simulating 
flow  in  large  3D  geometries  the  solver  has  been  parallelized  for  use  on  CPU  clusters 
using  MPI. 

2.3.3  Shock  Capturing 

One  of  the  greatest  restrictions  of  high-order  unstructured  solvers  is  their  inability  to 
handle  flow  discontinuities.  When  flows  develop  steep  gradients  such  as  shock  waves 
or  contact  surfaces,  non-physical  spurious  oscillations  arise  that  cause  the  simulations 
to  become  unstable.  For  higher-order  approximations,  it  is  typically  necessary  to  add 
explicit  dissipation  in  order  to  obtain  a  stable  solution.  But  this  has  a  negative  effect 
on  accuracy,  and  the  resolution  of  turbulent  scales.  The  development  of  numerical 
algorithms  that  capture  discontinuities  and  also  resolve  the  scales  of  turbulence  in 
compressible  turbulent  flows  remains  a  significant  challenge. 

In  this  study  the  shock  capturing  approach  of  Cook  and  Cabot  [21]  has  been  adapted 
to  the  SD  method.  The  scheme  adds  artificial  viscosity  to  the  flow  when  disconti¬ 
nuities  are  present  [11]  [12].  This  causes  the  discontinuities  to  be  smeared  out  over 
a  single  (high-order)  SD  element,  thus  precluding  the  development  of  any  spurious 
oscillations.  Results  illustrating  the  shock  pattern  for  supersonic  2D  flow  over  a  bump 
are  shown  in  Fig.  1. 
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(a) 


(b) 


(c) 

Figure  1:  Shock  pattern  resulting  from  2D  supersonic  flow  over  a  bump  calculated 
using  a  third-order  SD  method  with  artificial  viscosity.  The  computational  mesh  is 
shown  in  (a),  artificial  bulk  viscosity  is  shown  in  (b)  and  pressure  contours  are  shown 
in  (c).  Note  that  artificial  bulk  viscosity  is  only  added  in  the  presence  of  shocks. 
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2.3.4  Automatic  Mesh  Refinement 


An  optimal  computational  mesh  should  be  able  to  capture  details  of  the  flow  solu¬ 
tion,  yet  avoid  unnecessary  resolution  of  regions  in  which  the  flow  is  predominantly 
uniform.  In  order  to  generate  such  optimal  meshes  a  priori  knowledge  of  the  flow 
physics  is  required,  as  well  as  an  understanding  of  how  the  chosen  numerical  scheme 
performs.  Obtaining  such  information,  and  generating  such  optimal  meshes,  is  often 
very  difficult,  if  not  impossible.  Therefore,  to  mitigate  these  issues,  various  adaptive 
meshing  techniques  have  been  developed. 

Adaptive  meshing  techniques  allow  the  mesh  to  be  automatically  modified  based  on 
the  flow  solution.  The  nature  of  these  modifications  is  usually  determined  by  error 
indicators  calculated  from  the  flow  as  it  develops  (for  example  local  entropy  error 
may  be  used  as  an  indicator  for  isentropic  flows).  Regions  with  higher  error  are 
then  refined,  and  in  some  cases  regions  with  lower  errors  unrefined.  For  high-order 
schemes,  such  as  SD  and  DG  type  methods,  refinement  can  occur  in  a  variety  of 
ways.  These  include  reduction  of  the  element  size  ( /i-type  refinement),  increasing  the 
element  order  (p-type  refinement),  or  a  combination  of  both  (/ip-type  refinement). 

Various  error  estimation  and  mesh  refinement  strategies  have  been  implemented  and 
tested  within  our  in-house  2D  SD  flow  solver  [12],  An  illustration  of  how  /i-type  re¬ 
finement  has  be  used  to  improve  resolution  of  shock  waves  (forming  due  to  supersonic 
flow  over  a  bump)  is  shown  in  Fig.  2.  Note  that  artificial  viscosity  is  also  employed 
to  ensure  that  the  SD  scheme  is  stable  in  the  presence  of  the  shock  waves. 
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(a) 


Figure  2:  An  example  of  how  adaptive  h- type  mesh  refinement  can  facilitate  resolution 
of  shock  waves  at  Mach  1.4.  The  original  mesh  is  shown  in  (a)  and  the  original 
solution  in  (b).  The  automatically  refined  mesh  is  shown  in  (c)  and  the  resulting 
refined  solution  in  (d).  Note  that  artificial  viscosity  is  also  employed  to  ensure  that 
the  SD  scheme  is  stable  in  the  presence  of  the  shock  waves. 
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2.3.5  Mesh  Deformation 


Mesh  deformation  algorithms  have  been  implemented  within  our  in-house  2D  SD 
flow  solver.  In  the  schemes  we  have  implemented,  deformation  of  physical  bound¬ 
aries  causes  rigid  displacement  of  nearby  elements.  This  deformation  is  then  blended 
smoothly  into  the  mesh,  such  that  the  mesh  at  far  field  boundaries,  or  some  other 
desirable  portions  of  the  flow  domain,  remains  unaltered. 

Deformation  of  the  physical  mesh  is  achieved  via  time-dependent  variations  of  the 
metrics  and  the  Jacobian  that  define  the  mapping  of  each  physical  element  to  a 
standard  reference  element.  Such  an  approach  preserves  the  high-order  accuracy  of 
the  SD  method  since  the  governing  equations  are  always  solved  in  a  steady  reference 
element. 

An  example  of  mesh  deformation  and  mesh  blending  is  shown  in  Fig.  3.  The  undis¬ 
torted  passage  of  an  Euler  vortex  across  a  mesh  undergoing  significant  deformations 
in  shown  in  Fig.  4. 
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(a)  (b) 


Figure  3:  An  illustration  of  how  mesh  deformation  algorithms  can  be  used  to  blend 
a  rigid  displacement  of  elements  near  a  deformation  into  stationary  elements  at  the 
boundary  of  the  domain.  The  original  mesh  is  shown  in  (a)  and  the  resulting  deformed 
mesh  is  shown  in  (b). 


Figure  4:  The  initial  pressure  distribution  in  an  Euler  vortex  is  shown  in  (a).  The 
pressure  distribution  after  translation  over  a  mesh  undergoing  spatial  and  temporal 
deformations  is  shown  in  (b).  Note  the  the  structure  of  the  vortex  remains  uniform 
despite  the  non-uniform  mesh  deformation. 
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2.3.6  Convergence  Acceleration 


Efficient  p-rnultigrid  and  implicit  time  integration  methods  have  been  implemented 
to  accelerate  the  convergence  of  high-order  SD  schemes  to  steady-state  [14] .  To  test 
the  methods  2D  steady  viscous  flow  over  an  airfoil  at  Re  =  5000  has  been  investi¬ 
gated.  A  p-multigrid  approach  utilizing  a  3-level  V-cyclc  with  1-1-6- 1-1  smoothing 
iterations  was  compared  with  an  implicit  time-stepping  approach  utilizing  LU-SGS 
inner  iterations.  The  convergence  history  of  the  residual  is  shown  in  Fig.  5.  Results 
for  an  explicit  third-order  Runge-Kutta  scheme  are  also  shown  for  comparison.  It  can 
be  seen  that  the  p-multigrid  approach  reduced  the  computational  cost  by  a  factor  of 
eight  compared  with  the  explicit  third-order  Runge-Kutta  scheme,  while  the  LU-SGS 
approach  reduced  the  computational  cost  by  factor  of  100.  These  results  indicate 
that  by  using  efficient  convergence  acceleration  techniques,  the  computational  cost  to 
reach  a  steady-state  solution  using  the  SD  method  can  be  greatly  reduced. 
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Figure  5:  Plots  of  the  residual  against  CPU  time  for  an  SD  simulation  of  steady- 
state  viscous  flow  over  an  airfoil  at  Re  =  5000.  Results  for  three  time-integration 
schemes  are  shown.  These  include  a  p-multigrid  approach  utilizing  a  3-level  V-cycle 
with  1-1-6-1-1  smoothing  iterations,  an  implicit  time-stepping  approach  utilizing  LU- 
SGS  inner  iterations,  and  an  explicit  third-order  Runge-Kutta  scheme.  It  can  be 
seen  that  the  p-multigrid  approach  reduces  the  computational  cost  by  a  factor  of 
eight  compared  with  the  explicit  third-order  Runge-Kutta  scheme,  while  the  LU-SGS 
approach  reduces  the  computational  cost  by  factor  of  100. 


2.4  Code  Utilization 

2.4.1  Two-Dimensional  Flow  over  Cylinders 

Simulations  of  flow  over  a  pair  of  stationary  cylinders  have  been  undertaken  using  our 
in-house  2D  SD  viscous  compressible  flow  solver.  Vorticity  contours  for  cases  where 
the  cylinder  separation  is  three  times  the  cylinder  diameter  are  shown  for  various 
Reynolds  numbers  (Re)  in  Fig.  6.  For  further  details  see  Liang,  Premasuthan  and 
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Jameson  [15]. 
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Figure  6:  Plots  of  vorticity  over  a  pair  of  stationary  cylinders.  Results  are  presented 
for  Re  =  100  (a)  and  Re  =  200  (b). 

Simulations  of  flow  over  a  pair  of  counter  rotating  cylinders  have  also  been  undertaken 
using  our  in-house  2D  SD  viscous  compressible  flow  solver.  The  effects  of  Reynolds 
number,  compressibility,  and  rotation  speed  were  all  studied.  For  further  details  see 
Ou,  Liang,  Premasuthan  and  Jameson  [16]. 

2.4.2  Two-Dimensional  Plunging  and  Pitching  Airfoils 

Simulations  of  flow  over  plunging  and  pitching  NACA0012  airfoils  have  been  under¬ 
taken  using  our  in-house  2D  SD  viscous  compressible  flow  solver.  The  plunging  airfoil 
simulations  were  based  on  water  tunnel  experiments  performed  by  Jones  et  al.  [22] 
at  Re  =  1850.  Fig.  7  shows  that  the  pattern  of  vortical  structures  obtained  via  our 
numerical  simulations  compares  well  with  experimental  results.  In  particular,  the 
simulations  were  able  to  reproduce  the  fine  structures  occurring  in  the  wake  of  the 
plunging  airfoil. 
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Figure  7:  Vorticity  over  a  plunging  NACA0012  airfoil  at  Re  =  1850  calculated  using 
a  forth-order  SD  scheme  (a)  is  compared  with  an  analogous  experimental  result  from 
Jones  et  al.  [22]  (b).  The  numerical  results  compare  well  with  the  experimental 
results.  In  particular,  the  simulations  were  able  to  reproduce  the  fine  structures 
occurring  in  the  wake  of  the  plunging  airfoil. 


The  pitching  airfoil  simulations  were  based  on  water  tunnel  experiments  performed 
by  Koochesfahani  et  al.  [23].  Fig.  8  shows  two  comparisons  between  our  numerical 
simulations  and  experiments  (each  with  a  different  pitching  frequency).  In  the  first 
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case  the  wake  assumes  the  form  of  an  undulating  vortex  sheet,  and  in  the  second  case 


a  double-vortex  structure  is  seen  to  persist  downstream  of  the  airfoil.  In  both  cases 
results  of  our  simulations  are  seen  to  compare  well  with  the  experimental  data. 


Figure  8:  Vorticity  over  a  pitching  NACA0012  airfoil  at  Re  =  1.2  x  104  calculated 
using  a  forth-order  SD  scheme  (a)  is  compared  with  an  analogous  experimental  result 
from  Koochesfahani  et  al.  [23]  (b). 


2.4.3  Two-Dimensional  Interaction  of  a  Cylinder  Wake  with  a  Deforming 
Beam 

Simulations  of  flow  over  a  cylinder  connected  to  a  beam  undergoing  a  prescribed 
deformation  have  been  undertaken  using  our  in-house  2D  SD  viscous  compressible 
flow  solver.  The  setup  acts  as  a  simple  2D  model  of  a  body  connected  to  a  single 
flapping  wing.  Plots  of  vorticity  within  the  vicinity  of  the  configuration  are  shown  at 
various  instances  in  Fig.  9. 
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(a) 


(c) 


(d) 

Figure  9:  Plots  of  vorticity  over  a  cylinder  connected  to  a  beam  undergoing  a  pre¬ 
scribed  deformation.  Various  instances  in  time  are  shown;  specifically  the  peak  of  the 
upstroke  (a),  the  following  neutral  position  (b),  the  peak  of  the  downstroke  (c)  and 
the  following  neutral  position  (d). 
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2.4.4  Three-Dimensional  Turbulent  Channel  Flow 


Large  Eddy  Simulations  of  compressible  turbulent  channel  flow  were  undertaken  using 
the  SD  method  [17].  Contours  of  spanwise  vorticity  in  the  center  plane  of  the  channel 
obtained  using  a  fourth-order  SD  scheme  are  shown  in  Fig.  10.  The  predicted  mean 
and  root  mean  squared  velocity  profiles  were  found  to  be  in  good  agreement  with 
direct  numerical  simulation  results  obtained  by  Moser,  Kim  and  Mansour  [24], 


Figure  10:  Contours  of  spanwise  vorticity  in  the  center  plane  of  the  channel  obtained 
using  a  fourth-order  SD  scheme. 


2.4.5  Three-Dimensional  Turbulent  Flow  Over  an  Airfoil 

Simulations  of  flow  over  an  SD7003  airfoil  at  a  4°  angle  of  attack  have  been  performed 
using  our  in-house  3D  viscous  compressible  SD  flow  solver.  The  SD7003  airfoil  was 
selected  due  to  availability  of  existing  experimental  [25]  and  computational  data  [26] . 
Forth-order  accurate  simulations  were  undertaken  on  meshes  with  1.7  x  106  degrees 
of  freedom  at  Re  —  1  x  104  and  Re  =  6  x  104.  When  Re  —  lx  104  the  flow  was 
essentially  2D  with  close-to-periodic  vortex  shedding  (see  Fig.  11(a)).  The  computed 
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average  lift  coefficient  was  0.381  and  the  average  drag  coefficient  was  0.0497,  which 
are  consistent  with  those  obtained  by  Uranga  et  al.  [26].  At  a  Reynolds  number  of 
Re  =  6  x  104  transition  was  observed  to  take  place  across  a  laminar  separation  bubble 
(see  Fig.  11(b)). 

The  Q-criterion  Q  provides  a  mean  of  visualizing  vortex  cores  and  identify  turbulent 
structures.  It  can  be  calculated  as 


Q 


(flijflij  SijSij) 


(2.2) 


where 


--ij 


duj  _  duA 
dxj  dxi )  ’ 


^i\ 

dxi) 


(2.3) 


are  the  anti-symmetric  and  symmetric  parts  of  the  velocity  gradient  tensor  respec¬ 
tively.  Fig.  11  shows  instantaneous  iso-surfaces  of  Q  over  the  SD7003  airfoil  for  cases 
where  Re  =  1  x  104  and  Re  =  6  x  104. 
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(a) 


(b) 

Figure  11:  Instantaneous  iso-surfaces  of  Q  above  an  SD7003  airfoil  at  a  4°  angle  of 
attack.  When  Re  =  1  x  104  (a)  the  flow  is  essentially  2D  and  remains  laminar  over 
the  wing  surface,  while  at  Re  =  6  x  104  (b)  transition  takes  place  across  a  laminar 
separation  bubble. 


Use  of  the  SD  method  with  an  implicit  large  eddy  simulation  approach  appears  to  be 
capable  of  accurately  predicting  the  laminar  separation  and  transition  locations  over 
an  SD7003  airfoil  at  Re  =  6  x  104.  To  the  best  of  our  knowledge,  this  study  is  the 
first  attempt  to  analyze  transitional  flow  using  the  SD  method. 


2.4.6  Three-Dimensional  Flapping  Wing 

Preliminary  simulations  of  viscous  compressible  flow  over  a  SD7003  airfoil  undergoing 
a  prescribed  flapping  motion  have  been  undertaken  using  our  in-house  3D  viscous 
compressible  SD  flow  solver.  The  geometry  and  nature  of  the  prescribed  motion  are 
illustrated  in  Fig.  12. 
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Figure  12:  Deforming  SD7003  airfoil  at  bottom  of  cycle  (a),  middle  of  cycle  (b),  and 
top  of  cycle  (c). 


3  Conclusions 

A  summary  has  been  given  of  research  carried  out  under  award  number  FA9550- 
07-1-0195  from  the  Air  Force  Office  of  Scientific  Research.  The  objective  of  the  pro¬ 
posed  research  was  to  develop  and  implement  high-order  numerical  algorithms  for  the 
simulation  of  steady  and  unsteady  compressible  viscous  flows  with  shocks.  Notable 
achievements  have  been  made  in  three  areas,  namely  algorithm  analysis  and  develop¬ 
ment,  code  development,  and  code  utilization  (to  investigate  various  flow  problems). 
With  regards  to  algorithm  analysis  and  development,  it  has  been  proved  for  one¬ 
dimensional  linear  advection  that  the  SD  method  is  stable  for  all  orders  of  accuracy 
in  a  norm  of  Sobolev  type  (provided  that  the  interior  flux  collocation  points  are 
placed  at  the  zeros  of  the  corresponding  Legendre  polynomials).  Also,  a  new  range 
of  energy  stable  high-order  methods  based  on  the  so  called  FR  approach  have  been 
identified.  With  regards  to  code  development,  two-dimensional  and  three-dimensional 
compressible  viscous  flow  solvers  based  on  the  spectral  difference  method  have  been 


written.  Within  the  aforementioned  codes  a  range  of  shock  capturing,  automatic 
mesh  refinement,  mesh  deformation,  and  convergence  acceleration  algorithms  have 
been  implemented  and  tested.  With  regards  to  code  utilization,  viscous  compressible 
flow  over  various  2D  configurations  has  been  investigated.  These  configurations  in¬ 
clude  pairs  of  cylinders  (both  stationary  and  rotating),  plunging  and  pitching  airfoils, 
and  a  deforming  beam  in  the  wake  of  a  cylinder.  In  3D,  turbulent  channel  flow  and 
turbulent  flow  over  an  airfoil  have  been  investigated,  and  preliminary  simulations  of 
viscous  compressible  flow  over  a  flapping  wing  have  been  undertaken. 
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